#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_fit.h>
#include "cpgplot.h"

#define BUFSIZE 256

int main(int argc,char *argv[]){
  /* this code takes a spectral output from MIRIAD, that has
     frequency and amplitude, and computes the spectral index
     of that spectrum using the linear least-squares fitting
     routines from the GNU scientific library */
  FILE *spectrum=NULL;
  char line_in[BUFSIZE];
  double t_frequency,t_amplitude,c0,c1,cov00,cov11,cov01,sumsq;
  double *frequency=NULL,*amplitude=NULL;
  float *ffrequency=NULL,*famplitude=NULL,min_frequency,max_frequency;
  float min_amplitude,max_amplitude,*ffit=NULL;
  int n_frequencies=0,i;

  /* first we read in the file, which should be called spectrum.txt */
  spectrum=fopen("spectrum.txt","r");
  if (spectrum==NULL){
    printf("couldn't open spectrum.txt!\n");
    exit(-1);
  }
  min_frequency=1e6;
  max_frequency=-1e6;
  min_amplitude=1e6;
  max_amplitude=-1e6;
  while(fgets(line_in,BUFSIZE,spectrum)!=NULL){
    if (sscanf(line_in," %lf %lf",&t_frequency,&t_amplitude)==2){
      /* we read what we expected, add it to the array */
      n_frequencies++;
      frequency=realloc(frequency,n_frequencies*sizeof(double));
      amplitude=realloc(amplitude,n_frequencies*sizeof(double));
      ffrequency=realloc(ffrequency,n_frequencies*sizeof(float));
      famplitude=realloc(famplitude,n_frequencies*sizeof(float));
      frequency[n_frequencies-1]=log10(t_frequency);
      ffrequency[n_frequencies-1]=log10(t_frequency);
      if (frequency[n_frequencies-1]<min_frequency)
	min_frequency=frequency[n_frequencies-1];
      if (frequency[n_frequencies-1]>max_frequency)
	max_frequency=frequency[n_frequencies-1];
      amplitude[n_frequencies-1]=log10(t_amplitude);
      famplitude[n_frequencies-1]=log10(t_amplitude);
      if (amplitude[n_frequencies-1]<min_amplitude)
	min_amplitude=amplitude[n_frequencies-1];
      if (amplitude[n_frequencies-1]>max_amplitude)
	max_amplitude=amplitude[n_frequencies-1];
    }
  }
  fclose(spectrum);

  /* now print out the end points and the middle point */
  printf("frequency = %f amplitude = %f\n",frequency[0],amplitude[0]);
  printf("frequency = %f amplitude = %f\n",frequency[(int)(n_frequencies/2)-1],
	 amplitude[(int)(n_frequencies/2)-1]);
  printf("frequency = %f amplitude = %f\n",frequency[n_frequencies-1],
	 amplitude[n_frequencies-1]);

  /* and plot the spectrum */
  cpgopen("/xs");
/*   cpgenv(min_frequency,max_frequency,min_amplitude,max_amplitude,0,30); */
  cpgswin(min_frequency,max_frequency,min_amplitude,max_amplitude);
  cpgbox("BCLNTS",0,2,"BCLNTS",0,2);
  cpgline(n_frequencies,ffrequency,famplitude);
  cpglab("Frequency (GHz)","Amplitude (Jy)","");

  /* now do the fitting */
  gsl_fit_linear(frequency,1,amplitude,1,n_frequencies,&c0,&c1,&cov00,
		 &cov01,&cov11,&sumsq);
  ffit=malloc(n_frequencies*sizeof(float));
  for (i=0;i<n_frequencies;i++){
    ffit[i]=c1*ffrequency[i]+c0;
  }

  /* and print the result */
  printf("amplitude = %f x frequency + %f\n",c1,c0);
  /* and plot the result */
  cpgsci(3);
  cpgline(n_frequencies,ffrequency,ffit);
  cpgsci(1);

  cpgclos();
  exit(0);

}
